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1.  Introduction 


The  model  of  layered  random  medium  with  rough  interfaces  is  often  encountered 
in  many  applications  in  various  disciplines.  A  simple  approach  is  to  incoherently 
add  the  contributions  of  volumetric  and  surface  fluctuations  [36,  16].  However,  this 
is  valid  only  when  we  are  in  the  single  scattering  regime  [7,  20].  There  are  some 
other  hybrid  approaches  [23,  6]  which  take  into  consideration  some  multiple  scat¬ 
tering  effects.  Brown  [3]  outlines  an  iterative  procedure  which  properly  includes  all 
multiple  scattering  interactions.  However,  it  does  not  appear  feasible  to  carry  out 
the  calculation  beyond  one  or  two  iterations.  Among  the  other  methods  currently 
used,  perhaps  the  most  widely  used  approach  is  the  radiative  transfer  (RT)  approach 
[31,  14,  13,  27,  4,  32,  17,  10].  Here  one  formulates  the  scattering  and  propagation  in 
each  layer  by  using  the  radiative  transfer  equation  which  involves  only  the  parameters 
of  the  medium  of  that  layer.  The  boundary  conditions  are  derived  separately  and  in¬ 
dependently  using  some  asymptotic  procedure  developed  in  rough  surface  scattering 
theory  [2,  1,  33].  The  RT  equations,  along  with  the  boundary  conditions,  comprise 
the  system  that  describes  the  problem. 

Although  this  procedure  appears  to  be  reasonable  and  sound  it  is  apparent  that 
certain  approximations  are  involved  and  we  would  like  to  know  the  conditions  under 
which  this  kind  of  approach  is  appropriate.  One  way  to  better  understand  the  RT 
approach  is  to  compare  it  with  the  more  rigorous  wave  approach.  For  the  case  of 
unbounded  random  media  it  was  found  [35]  that  the  RT  approach  is  applicable  when 
(a)  one  uses  the  quasi-uniform  field  approximation,  (b)  one  uses  the  ladder  approxi¬ 
mation  to  the  intensity  operator  of  the  Bethe-Salpeter  equation,  and  (c)  the  medium 
is  statistically  quasi-homogeneous. 

However,  our  problem  has  bounded  structures  which  are  randomly  rough.  There¬ 
fore  it  remains  to  be  seen  whether  the  conditions  arrived  at  in  the  case  of  unbounded 

random  media  will  be  sufficient  for  our  problem. 
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In  our  project  we  employed  a  statistical  wave  approach  using  surface  scattering 
operators  [33,  21]  to  derive  the  transport  equations  for  our  multilayer  problem.  In 
this  process  we  found  that  there  are  more  conditions  implied  when  we  choose  to  apply 
the  RT  approach  to  our  problem  than  it  is  widely  believed  to  be  necessary.  One  such 
condition  is  the  weak  surface  correlation  approximation.  This  means  that  the  RT 
approach  places  certain  restrictions  on  the  type  of  rough  interfaces  that  it  can  model 
accurately. 

This  report  is  organized  as  follows.  In  Section  II  we  describe  the  geometry  of  the 
problem.  In  the  next  section  we  describe  the  RT  approach  to  the  problem.  In  Section 
IV  we  describe  the  wave  approach  to  this  problem.  In  Section  V  we  transition  to  the 
transport  equation  system.  The  report  concludes  with  a  discussion  of  our  findings. 

2.  Description  of  the  Problem 

The  geometry  of  the  problem  is  shown  in  Figure  1.  We  have  an  TV-layer  random 
medium  stack  with  rough  interfaces.  The  permittivity  of  the  j-th  layer  is  6j  +  lj  (r) 
where  6j  is  the  deterministic  part  and  £j  is  the  randomly  fluctuating  part.  The 
permeability  of  all  the  layers  is  that  of  free  space.  The  randomly  rough  interfaces 
are  given  as  z  =  Zj  +  Cj(r_i_)-  R  is  assumed  that  £j  and  Q  are  zero-mean  isotropic 
stationary  random  processes  independent  of  each  other.  Thus,  on  the  average  the 
interfaces  are  parallel  planes.  Let  zq  =  0,  and  let  dj  be  the  thickness  of  the  j-th  layer. 
The  media  above  and  below  the  stack  are  homogeneous  with  parameters  eo,  &o,  and 
ejv+i,  i,  respectively.  This  system  is  excited  by  a  monochromatic  electromagnetic 
plane  wave  and  we  are  interested  in  formulating  the  resulting  multiple  scattering 
process. 

3.  Radiative  Transfer  Approach 

Multiple  scattering  in  a  complex  environment  is  well  described  by  the  radiative 

transfer  theory.  This  theory  is  not  only  conceptually  simple  but  also  very  efficient. 
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Figure  1:  Geometry  of  the  problem. 


The  fundamental  quantity  here  is  the  specific  intensity  I  which  is  governed  by  the 
following  equation  [5,  28,  12] 


(1) 


One  may  regard  this  equation  as  a  statement  of  conservation  of  energy  density  I 
which  is  a  phase-space  quantity  at  position  r  and  direction  s.  7  is  the  extinction 
matrix  which  is  a  measure  of  loss  of  energy  in  direction  s  due  to  scattering  in  other 
directions.  P  is  the  phase  matrix  representing  increase  in  energy  density  in  direction 
s  due  to  scattering  from  neighbouring  elements,  is  the  solid  angle  subtended  by  s. 
Given  the  statistical  characteristics  of  the  medium  one  can  readily  calculate  the  phase 
matrix.  The  extinction  matrix  is  hence  calculated  using  the  relation  7=  j  P(s',  s)dfl'. 
The  specific  intensity  in  each  layer  is  governed  by  an  equation  similar  to  (1).  Since 
our  layer  problem  has  translational  invariance  in  azimuth  the  RT  equation  for  the 
m-th  layer  takes  the  following  form, 


COS  0  ^  ^  Im  (.2,  s)T  7  m  I  m  dll  Pm(®i  ■S  )I  m  (%  j  S  ) 

m 
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where  the  subscript  m  denotes  that  the  quantity  corresponds  to  those  of  the  m-th  layer 
and  d  is  the  elevation  angle  of  s.  This  set  of  RT  equations  is  complemented  by  a  set  of 
boundary  conditions  which  are  in  turn  based  on  energy  conservation  considerations. 
To  be  more  precise,  we  impose  the  condition  that  the  energy  flux  density  at  each 
interface  is  conserved.  This  leads  to  the  following  boundary  conditions  on  the  m-th 
interface 


rm(zm,  s)  =  f  <m'  <ftm+1,m(s,  s'))  I dm(zm,  s')  +  J  dD!  (Tm,m+1(s,  s'))  1“  +1(zm,  s').  (3) 

The  boundary  conditions  on  the  (m  —  l)-th  interface  are  given  as 


ldm{zm-i,s)  =  /  (IQ'  (Km_1}rn(s,s'))Ill(zm_1,s')+  /  dO'  (Tmm_i(s,  s'))  s'), 


(4) 


where  TZmn  and  Tmn  are  the  local  reflection  and  transmission  Miiller  matrices.  To 
be  more  specific,  TZmn  represents  the  reflection  Muller  matrix  of  waves  incident  from 
medium  n  on  the  interface  that  separates  medium  m  and  medium  n.  The  superscripts 
u  and  d  indicate  whether  the  intensity  corresponds  to  a  wave  travelling  upwards  or 
downwards.  The  integrations  in  these  expressions  are  over  a  solid  angle  (hemisphere) 
corresponding  to  s'.  Suppose  we  have  a  plane  wave  incident  on  this  stack  from  above. 
Then  the  downward  travelling  intensity  in  Region  0  is 


Io(z,s)  =  B 05  (cos  d0  -  cos  d*)  d(0o  -  00, 


(5) 


where  B0  is  the  intensity  of  the  incident  plane  wave  and  {d,,  0*}  describes  its  direction. 
Since  there  is  no  source  or  scatterer  in  Region  N  +  1, 


=  °- 

Notice  again  that  these  boundary  conditions  represent  conservation  of  intensity  at 
the  interfaces. 

We  should  point  out  that  the  radiative  transfer  approach  as  applied  to  a  particular 

problem  is  only  a  model  based  on  certain  assumptions.  Since  the  RT  theory  is  used  in 
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a  variety  of  applications,  the  particular  assumptions  involved  are  described  in  terms 
of  different  terminologies,  specific  to  the  discipline  where  it  is  used.  One  good  way  to 
understand  in  more  general  terms  the  RT  approach  and  the  underlying  assumptions 
is  to  compare  it  with  the  more  rigorous  statistical  wave  approach.  For  the  case  of  an 
unbounded  random  medium  this  kind  of  study  was  carried  out  in  the  seventies  [35]. 
From  that  study  we  learn  that  the  radiative  transfer  theory  can  be  applied  under  the 
following  conditions: 

1.  Quasi-stationary  field  approximation. 

2.  Weak  fluctuations: 

(a)  Ladder  approximation  to  the  intensity  operator. 

(b)  Kraichnan  approximation  to  the  mass  operator. 

3.  Statistical  homogeneity  of  the  medium  fluctuations. 

These  are  the  well-known  conditions  that  we  associate  with  the  RT  approach. 
However,  our  problem  has  bounded  structures  and,  further,  they  are  randomly  rough. 
The  question  is  this:  are  the  above  conditions  sufficient  to  apply  the  RT  approach  for 
our  problem?  This  is  the  motivation  for  our  investigation.  We  follow  the  statistical 
wave  approach  to  this  problem,  derive  the  equations  for  the  coherence  functions,  and 
hence  make  the  transition  to  the  RT  equations.  This  procedure  enables  us  to  better 
understand  the  necessary  conditions  for  using  the  RT  approach  for  our  problem. 

4.  Wave  Approach 

The  following  are  the  equations  that  govern  the  waves  in  the  layer  structure: 

V  x  V  x  E j  —  kjEj  =  VjEj  j  —  1,  ■  ■  • ,  N,  (6) 

where  Vj  =  r)  represents  the  volumetric  fluctuation  in  Region  j.  For  the 

homogeneous  regions  above  and  below  we  have 

VxVxEq  -  /cqE  o  =  0,  (7a) 
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(7b) 


V  x  V  x  Ej\r+i  -  fcjv+i  Ejv+i  —  0. 
The  boundary  conditions  at  the  j-th  interface  are 


h  x  Ej(rj_,  0)  =  n  x  EJ+i(r±,  G), 


(8a) 


n  x  V  x  E^rj^Ci)  =  n  x  V  x  Ei+i(r_L,  (,•)■ 


(86) 


where  h  is  the  unit  vector  normal  to  the  j-th  interface  with  normal  pointing  into  the 
medium  j .  This  system  is  complemented  by  the  radiation  conditions  well  away  from 
the  stack.  We  assume  that  we  know  the  solution  to  the  problem  without  volumetric 
fluctuations,  and  denote  it  as  E.  The  corresponding  Green’s  functions  denoted  as  G ^ 
are  governed  by  the  following  set  of  equations: 


V  x  V  x  Gjk(r,  r')  -  kjGjk(r,  r')  =  I5jk5(r  -  r),  (9a) 

n  x  Gjk(r±,  Cf,  r')  =  nx  G{j+1)k(r±,  Q;  r),  (9b) 

n  x  V  x  Gjk(r±,  Q]  r)  =  n  x  V  x  G(i+i)fc(rj_,  Q]  r')-  (9c) 


Another  pair  of  equations  similar  to  (9b)  and  (9c)  corresponding  to  (j  —  l)-th  interface 
need  to  be  added  to  this  list.  I  here  is  unit  dyad.  Using  these  Green’s  functions  and 
the  radiation  conditions  the  wave  functions  can  be  represented  as 

N  r 

Ej(r)  =  Ej(r)  +  ^  /  dr'Gjk(r,r')vk(r')Ek(r')  j  =  0, 1,  •  •  • ,  N  +  1  (10) 

fe=1 

where  Qk  =  {r'±;  (k  <  z'  <  Cfc-i}-  Note  that  v0  =  vn+i  =  0.  We  first  average  (10) 
w.r.t.  volumetric  fluctuations  to  get 
<E,(r)).  =  E,(r)+ 

*'  /  dr"Gjk( r1r1(Gll(r',0),h(r>M)(El(r")), 

•k 

(11) 


where  G ki  is  governed  by  the  following  system  of  equations 


V  x  V  x  Gki(r,  r)  -  hkGki(r,  r')  =  lSki§(r  -  r)  +  vkGH(r,  r') 
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n  x  Gfc/(r_j_,  (k;  r)  =  n  x  G{k+1)i(r±1  (k,  r), 

n  x  V  x  Gw(r_L,  (k,  r')  =  n  x  V  x  G(fc+i)Z(r_L,  (k;  r'). 

Here  h  being  the  unit  vector  normal  to  the  fc-th  interface.  We  also  have  a  similar  set 
of  boundary  conditions  on  the  k  —  1-th  interface.  The  subscript  v  is  used  to  denote 
averaging  with  respect  to  volumetric  fluctuations.  Here  we  have  used  a  first  order 
approximation  to  the  mass  operator  based  on  weak  fluctuations.  The  volumetric 
fluctuations  in  different  regions  are  assumed  to  be  uncorrelated,  which  means  that 

{vk(r')vi( r"))  =  8klCk(r'  -  r"),  (12) 


where  Ck  is  the  correlation  function  of  the  volumetric  fluctuations  in  Region  k.  We 
have  assumed  that  the  fluctuations  of  the  parameters  of  our  problem  are  Gaussian 
and  statistically  homogeneous.  Inserting  (12)  in  (11)  and  employing  V  x  V  x  I  —  k? 
on  (11)  we  get 

V  x  V  x  (EJ(r))„  -  kj(Ej(r))v  =  I  dr'( G^r,  r')>„C,-(r  -  r,)(Ej(r'))..  (13) 

Next  we  average  (13)  over  the  surface  fluctuations, 


V  X  V  X  (E,(r)>„,  -  k]{ E,(r)) 


*'(<Gjj(r.r'))»cj(r  “  r'){Ej(r'))„)^ 


where  the  subscript  s  denotes  averaging  over  surface  fluctuations  and  Qj 
z’  <  Zj_i}.  We  approximate  ( (Gjj(r, r')),,^^  —  r/)(EJ(r/))^\  as 


\  IS 

(Gjj(r,r'))vsCj(r  -  r')(Ej(r/))va  and  obtain 


(r±;%  < 


V  x  V  x  (Ej(r))„-k?(Ej(r))„  =  [  dr'( Gfl(r,r'))MCj(r-r')(Ej(r'))M.  (14) 

We  call  this  the  weak  surface  correlation  approximation.  We  will  later  find  that  this 

is  one  additional  approximation  necessary  to  arrive  at  the  RT  system.  First  note  from 

(14)  that  (V  x  Vx  I  —kj){Ej(r))va  =  0  for  j  =  0,  N+ 1.  This  means  that  the  coherent 

propagation  constants  in  regions  above  and  below  the  layer  stack  are  unaffected  by 

the  fluctuations  of  the  problem.  However,  they  indeed  get  modified  within  the  stack 
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region.  On  writing  (14)  as  (V  x  V x  I  —k2  —  £)(^-)  =  0,  where  C  denotes  the  integral 
operator  dr'  ( Gjj(r ,  r'))vs  Cj(r  —  r'),  we  infer  that  Xj  =  Vitf  +  £)  represents  the 
mean  propagation  constant  in  f Observe  that  Xj  depends  explicitly  on  the  volu¬ 
metric  fluctuations  in  Region  j  and  implicitly  on  the  fluctuations  of  the  stack,  both 
volumetric  and  surface.  This  is  in  contrast  to  the  RT  approach  where  7 j  depends 
exclusively  on  the  volumetric  fluctuations  in  Region  j.  Moreover,  Xj  depends  on  the 
polarization  if  the  fluctuations  of  the  problem  are  anisotropic.  Further,  even  if  the  vol¬ 
umetric  fluctuations  are  isotropic  Xj  will  be  polarization-dependent  because  of  surface 
reflections.  This  is  in  contrast  to  the  RT  approach  where  7 j  is  polarization-dependent 
only  when  the  volumetric  fluctuations  are  anisotropic.  Therefore  the  question  is  this: 
when  do  the  effects  of  boundaries  on  the  mean  propagation  constants  become  negli¬ 
gible?  A  first  order  solution  to  the  above  dispersion  relation  shows  that  in  situations 
where  the  thickness  of  the  layer  is  larger  than  the  corresponding  mean  free  path  the 
influence  of  the  boundaries  on  the  mean  propagation  constants  become  negligible,  as 
in  the  case  of  the  RT  system. 

Since  the  problem  is  invariant  under  translations  in  azimuth  the  mean  wave  func¬ 
tions  for  our  problem  have  the  following  form: 

(E)'(r))„s  =  exp(ik_u  ■  r){A^(k_u)pt  exp [itfjz]  +  B^(k±i)pj  exp 

j  =  1,2, •••,7V,  (15) 

(Eo(r))„s  =  exp(7kj_j  •  r)  jp^  exp  [~ik0ziz\  +  ^(k^p^  exp[ik0ziz\j ,  (16) 

and 

(E^+1(r))„s  =  exp(7kj_j  •  r)Tp(k_Li)p~+1  exp[-ik{N+1)ziz\,  (17) 

where  the  superscript  p  stands  for  the  polarization,  either  horizontal  or  vertical,  p  is 
the  unit  vector  representing  polarization,  qj  is  the  ^-component  of  Xj-  The  subscript 
i  is  used  to  indicate  that  the  wave  vector  is  in  the  incident  direction.  R  and  T  denote, 
respectively,  the  mean  reflection  and  transmission  coefficient  of  the  stack.  Aj  and  Bj 
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denote,  respectively,  the  mean  coefficients  of  up-going  and  down-going  waves  in  the 
j'-tli  layer.  Based  on  this  we  can  formulate  the  waves  averaged  w.r.t.  volumetric 
fluctuations  as 


<E, 


(r))p  =  —  / 

1  )/v  4?r2  J 


dk±  exp(ikj_-r)  j.AP9(kj_,  k_u)q+  exp[iqjz]+B^q(k_L,  ku)qj  exp [-iqjz 

j  =  (18) 


(E0(r))p  =  exp(ik_u  •  r)  exp[-ik0ziz]po  + 

+  4^2  /  ^k-LexP(^k_L  •  r)i?P9(k_L,k_Li)q(]'exp[iA;o22:],  (19) 

and 

(EwW)^^  J  dk±  exp(fk^  •  r)TP9(k^,  k_u)qjv+i  exp[— ik(N+i)zz\.  (20) 

where  Aj,  Bj,  R  and  T  are  now  integral  operators  representing  scattering  from  rough 
interfaces.  The  boundary  conditions  associated  with  the  above  equations  at  the  j-th 
interface  are 


“x(Ei(ri,(i))„  =  nx(Ei+i(ri,(i))#  J  =  1,2,---,  iV  (21a) 

and 


n  x  V  x  (Ej(ri,Ci))„  =  n  x  V  x  (EJ+1(r±Xj))v  j  =  1,  2,  ■  ■  ■ ,  N.  (21  b) 

The  above  system  may  be  solved  either  numerically  or  by  anyone  of  the  asymptotic 
methods*  available  in  rough  surface  scattering  theory  [2,  1,  33]  to  evaluate  the  mean 
coefficients  that  appear  in  (15)-(17). 

We  proceed  now  to  the  analysis  of  the  second  moments,  by  starting  with  (10).  For 
convenience  we  write  it  in  symbolic  form  as 

N 

E j  =  E j  +  Gjkvk Efc.  (22) 

k= i 

We  take  the  tensor  product  of  this  equation  with  its  complex  conjugate  and  average 
w.r.t.  volumetric  fluctuations  and  obtain 

N  N  N  N 

(Ej  ®  E*)^  =  (E,-),,® (E*)v  +  J2  (Et  (8)  E*,)v  ,  (23) 

k=  1  k'= 1  1=1  l'= 1 

*It  is  necessary,  however,  to  meet  the  weak  surface  correlation  approximation  employed  earlier 
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where  K  is  the  intensity  operator  of  the  volumetric  fluctuations.  Employing  the  weak 
fluctuation  approximation  we  approximate  K  by  its  leading  term 

K kk'iv  —  {vk  ®  vl)Skk'u /  I  .  (24) 

Next  we  average  (23)  w.r.t.  the  surface  fluctuations  and  employ  the  weak  surface 
correlation  approximation,  as  before,  to  get 

N 

<E;  «  E,*>„  =  ((E,)„  ®  <E*}„)  +  Y.  (  (&,*>„  ®  (G**>„  )  (»k  ®  Vk)  { E*  ®  E;>„  . 

S  k= 1 

(25) 

The  above  is  an  equation  for  the  second  moment  of  the  wave  function  E,  which  can 
be  decomposed  into  a  coherent  part  E  and  a  diffuse  part  E.  Therefore, 

(E®E*)  =  (E)  <g>  <E*)  +  (E0E*).  (26) 

The  coherent  part  is  not  of  much  interest;  we  know  that  it  is  specular  for  our  problem. 
The  diffuse  or  the  incoherent  part  is  of  more  interest.  Therefore  we  write  (25)  in  terms 
of  diffuse  fields: 

_ _  _ _  N 

(e,®E-)  =  ((Ej)„®{E*)v)  +5]  (27) 

S  k= 1 

where  (Ej)^  is  the  fluctuating  part  of  (Ej)^.  Let  us  now  write  (27)  in  more  detail  as: 

at 

(E,(r)  ®  E*(r'))  =  ((E,(r))„  ®  (E‘(r'))„)  +£  dr,  dr[ 

S  k=l  J 

(  (Gjk(r,ri))v  ®  <G*fc(r',r'1)>^jA;fc|4C'fe(r1  -  r'J)  (Efc(ri)  ®  E*^))^  .  (28) 

As  it  stands  this  equation  is  not  convenient  for  seeking  solution  either  analytically  or 

numerically.  Besides,  one  important  goal  for  us  is  to  investigate  the  conditions  needed 

for  employing  the  radiative  transfer  approach.  With  this  in  mind  we  introduce  Wigner 

transforms.  Note  that  (28)  is  an  equation  for  the  coherence  function.  On  the  other 

hand  the  RT  equation,  as  we  saw  earlier,  is  an  equation  for  the  specific  intensity 

which  is  a  ‘phase-space’  quantity.  Wigner  transforms  serve  as  a  bridge  to  link  these 

two  quantities  [34,  9,  18,  24], 
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We  introduce  Wigner  transforms  of  waves  and  Green’s  functions  as 


r  +  r' 


,k 


£ 

L'n 


ri  +  r1 


(U^.k)  =  J d(r  —  r')  (Em(r)  ®  E^(r'))  e-'k'<r_r')> 

■  l)  =  j  i(r-r')  J rfi'r  -  r',  j 

g-ilc  fr-r'Jgil  frj-ri,/  (Gmn(r,ri)>„  ®  (6^,^,^)), 


(29) 


(30) 


In  terms  of  these  transforms  (28)  becomes 

N 


£n(r,k)  =  ^(r,k)  +  -^^^|/cnjl  £  dr'  I  da  j 


11  1  dr’  I  da  I  d(3  Gmn(r,  k|r',  a)$n(a-f3)£n(r' ,  (3), 

(31) 


where  is  the  spectral  density  of  the  volumetric  fluctuations  in  the  n-th  layer.  We 
have  used  the  superscript  s  in  the  first  term  to  indicate  that  this  is  due  to  surface 
scattering  as  defined  by  the  first  term  in  (28). 

The  fact  that  our  problem  has  translational  invariance  in  azimuth  implies  the 
following: 


£m(r,k)  =  £m(z,  k), 

Gmn( r,  k|r',  1)  =  gmn(z ,  k\z',  1;  r±  -  r^). 


(32a) 

(326) 


Using  these  relations  in  (31)  we  have 

£m(z,  k)  =  £‘Jz,  k)  +  £  \kn\‘  J°"~‘  dz’  J  da  J  d/3 

Gmn(z ,  k| z\  a ;  0)<Fn(a:  -  f3)£n(z',  f3),  (33) 

where  Gmn{z,  k|U,  ck;  0)  is  the  Fourier  transform  of  Gmn{z,  k|U,  ck;  rj_  —  rj_)  w.r.t.  r  !.- 
r'j_  evaluated  at  the  origin  of  the  spectral  space.  To  proceed  further  we  need  to 
evaluate  Qmn-  Furthermore,  we  need  to  relate  this  system  with  that  of  RT,  which 
involves  the  boundary  conditions  at  the  interfaces.  In  view  of  this  we  need  to  identify 
the  coherence  functions  corresponding  to  up-  and  down-going  wave  functions.  To 
facilitate  this  we  decompose  (Gmn)v  into  its  components, 


(Gxmn/v  ^mn^m 


■  puu  _i_  r-^ud 
'  ^ mn  '  ran 


Gdu  |  p 
m.n.  ~T 


dd 

mn  i 
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(34) 


where  the  first  term  is  the  singular  part  of  the  Green’s  function.  The  superscripts  u 
and  d  indicate  up-  and  down-going  elements  of  the  waves.  The  other  components  are 
due  to  reflections  from  boundaries.  These  are  formally  constructed  using  the  concept 
of  surface  scattering  operators  as  follows  [33], 

<Gl(r,r')>r  = 


1 


(2tt) 


dk±  /  dkl{^n(k±,k/±)}^eiki'r+“^(kj-)2e-^'r'-^(k^)2',  (35) 


where  S^n  is  the  surface  scattering  operator.  The  superscripts  a  and  b  on  S  are 
used  to  indicate  whether  the  waves  are  up-going  or  down-going.  In  the  exponents, 
a,b  =  1  if  the  waves  are  up-going.  We  let  a,  b  =  —  1  if  the  waves  are  down-going.  The 
^-component  of  the  mean  propagation  constants  in  the  n-th  layer  is  denoted  as  qn. 

We  recall  that  Qmn  is  the  Wigner  transform  of  ^  ( Gmn)v  ®  ^  .  The  super¬ 

scripts  /i,  v  stand  for  polarization,  either  h  or  v.  It  is  important  to  note  that  only  the 
in  the  quasi-uniform  limit  does  the  Wigner  transform  of  the  coherence  function  lead 
to  the  specific  intensity  of  the  RT  equation.  For  our  layer  geometry  the  Green’s  func¬ 
tion  is  nonuniform.  However,  each  of  its  components  given  in  (34)  is  quasi-uniform. 
When  we  use  (34)  to  perform  the  Wigner  transform  we  ignore  all  cross  terms.  In 
other  words,  we  make  the  following  approximation, 


C  ~  A  C°  -I-  C]uu  I  r*ud  _i_  r>du  r>dd 

ymn  —  (Jmnyrn  '  ^ran  '  ^ran  ' 

where  Q^n  is  the  Wigner  transform  of  ^  ^  •  Most  of  the  cross  terms 

are  nonuniform  and  may  be  neglected  under  the  quasi-uniform  field  assumption.  A 
couple  of  cross  terms  are  quasi-uniform  and  their  inclusion  lead  to  phase  matrices 
that  are  different  from  those  of  the  RT  system.  It  turns  out  that  such  additional 
coherence  terms  become  negligible  when  the  layer  thickness  is  of  the  same  order 
or  greater  than  the  mean  free  path  of  the  corresponding  layer.  It  is  under  these 
conditions  the  approximate  expression  for  Qmn  given  above  is  good. 

With  the  introduction  of  this  representation  for  Qmn  in  (31)  we  can  trace  up-  and 
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down-going  waves  to  obtain  the  following  equations  for  the  coherence  function: 


£“  fok)  =£“(*,  k)+ 

1 


f  ^|fcm|4  J  dz'  j  da  I  dp  £>(;?,  k|z' ,  a;0)$m(a-  P)Sm(z' ,(3) 


N 


(2ir) 


£ 1*” 

71=1 


f*zn  —  1 


£ttz,k)  =  e£{z,k)+ 


f*zm—  1 


(2tt) 


\kr. 


dz'  da  d/3  Q^n(z,  k|z',  a;  0)<f>n(a  —  P)£“(z',  (3) 

Zn  J  J 

(36  a) 

dz'  [  da  [dp  Gm(z,  k|^/,a;0)$TO(a  -  f3)£m(z',/3) 


N 


(2vr)6  ^ 

v  '  n= 1 


V  |fcn 


dz'y  da  y  dp  Q^n(z,\<L\z' ,a\0)(&n(a  -  P)£“(z',P) 

(36 6) 


Note  that  summation  over  a  =  {w,  d}  is  implied  in  the  above  equations.  The  first 
term  in  these  equations,  Ssa,  represents  the  contribution  due  exclusively  to  surface 
scattering,  and  has  the  following  form: 


k)}'  =  2tt<S  <j  kz  -  K(k±)  +  C(k±)] 


0ia\q™-<im}z  /  i  Y/Ja  ^Y!,a*''Ul/ 


(k±,k±i)>  (37) 

/  S 

where  £^  is  the  amplitude  of  the  up-going  wave  in  the  m-th  layer  after  volumetric  av¬ 
eraging  is  performed.  This  means  that  it  is  a  random  function  of  surface  fluctuations. 
When  we  substitute  (37)  and  the  expressions  for  Qrnn  in  (36)  we  find  that 


Cfe k) }  =2*l\k,-  bl*( ki)  +  C(kJ]  \  {  £“(z, kj 


(38) 


/ IV  I  2  J  k  J  [1U 

On  substituting  this  in  (36)  and  differentiating  w.r.t.  z  we  obtain  the  following 
transport  equations: 


^  -  *  [?/x(k±)  -  q*v{ k±)]  |  £™u(z,  kj_)  =  Slv(z,  kj_)  +  ^ 


(2 WJ  d<X±Sl >X* 


<j  k±  -  a±;  ^  [gM(kj_)  +  q*( kj_)]  -  \  [q^>{a±)  +  <?*/(«l)]  }  {fi-fi,){u-u')E^vl(z,  a±), 


(39a) 
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J:  -  i  Mk±)  -  Qt( k±)] }  k±)  =  k±)  +  l^ja  /  da±S»Sv* 

®m  |k±  -  a±;  Mk±)  +  g*(k±)]  -  ^a[gM'(ajJ  +  g*,(a±)]J  {n-^){iy-u')£^v,(z,  a±), 

(396) 

~  a 

where  £  represents  scattering  due  to  the  coherent  part  of  £,  whereas  the  integral 
term  in  (39)  represents  scattering  due  to  the  diffuse  part  of  £.  We  may  also  regard 

—CL 

£fll/  as  the  source  to  our  transport  equations;  it  is  given  as 

=  \km\  !4,„  jkj.  -  k±1:  \  [qljk_]  -  ri'jkj)  -  h>  [j„(k±j)  +  ?;(k±i)|  1 

(s“S")  (m  ■  M.X"  '  Vi)  ((SL)m^)  ({S^UE«y.  (40) 

where  summation  over  6  is  implied.  Note  that  £a  in  (39)  includes  both  £u  and  £d 
(corresponding  to  up-  and  down-going  waves).  When  the  superscripts  {a, 6}  corre¬ 
spond  to  u,  the  value  of  {a,  6}  in  the  argument  of  <3>m  takes  the  value  +1;  on  the  other 
hand,  when  the  superscripts  {a,  6}  correspond  to  d  the  value  of  a  in  the  argument  of 
takes  the  value  —1.  Since  all  quantities  in  (39)  and  (40)  correspond  to  the  same 
layer  m  we  have  dropped  the  subscript  m  to  avoid  cumbersome  notations.  Summa¬ 
tion  over  [i!  and  v'  is  implicit  in  (39).  Similarly,  summation  over  p*  and  zy  is  implicit 
in  (40).  Sf^Q  is  the  scattering  amplitude  of  waves  with  direction  6  in  m-th  layer  due 
to  wave  incident  in  Region  0.  To  obtain  appropriate  boundary  conditions  we  have 
to  go  back  to  the  integral  equation  representations  for  £™u  and  £dv  and  observe  their 
behaviour  at  the  interfaces  and  try  to  find  a  relation  between  them. 

After  some  manipulations  we  arrived  at  the  following  boundary  conditions.  At  the 
(m  —  l)-th  interface  we  have 

£m(Zm-l,  kj_)  =  J  dk'±  (iZm-i,m  (kj.,  k'±)^  £“  (zm_i,  k'j.),  (41a) 

with  7Z=IC  <S>  R  where  R,n_ijm  is  the  stack  reflection  matrix  (not  the  local  reflection 
matrix)  for  a  wave  incident  from  below  on  the  (m  —  l)-th  interface.  Similarly 

i,kj_)  =  J  dk'±  (nm+i,m  (k±,  k'±) ^  £* (zm_i ,  k'±) ,  (416) 
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where  lZm+i,m  is  the  tensor  product  of  stack  reflection  matrix  for  a  wave  incident 
from  above  on  the  m-th  interface. 

We  were  able  to  obtain  the  boundary  conditions  only  after  imposing  certain  ap¬ 
proximations  as  given  below.  Consider  the  following  identity: 

Si,  =  F„,  {S>  +  Si,}  F„,  (42) 

where  Fm  =  diag  eiqvdm  }.  Notice  that  this  is  an  operator  relation  where  all 

elements  are  operators.  Taking  the  tensor  product  of  (42)  with  its  complex  conjugate 
we  have 

si,®si  =  (f,„  0  f:„)  ®  fC_i,m)  ({s’ + s~m}  ®  {§> + sr„n  (f„  ®  f^ 

(43) 

Next  we  average  (43)  w.r.t.  surface  fluctuations  and  get 

(si,®s** )  ~  (F,„  0  F;„)  (fL_1jTO  ®  ( {§i  +  S“„,}®{s’  +  Si,}* )  (Fm  ® 

(44a) 

where  we  the  two  tensor  products  in  the  middle  are  assumed  to  be  weakly  correlated. 

A  further  approximation  that  we  impose  is  given  as  follows 

( {S’  +  Si,}  ®  {S>  +  Si,}* )  ~  S’  0  S’*  +  (Si,  ®  Si).  (446) 

This  is  similar  to  the  approximation  we  used  while  computing  the  Wigner  transforms 
of  the  Green’s  functions.  These  are  the  kinds  of  approximations  required  to  arrive  at 
our  boundary  conditions. 

5.  Transition  to  Radiative  Transfer 

Now  we  have  to  transition  from  this  transport  equation  (44)  to  the  phenomeno¬ 
logical  radiative  transfer  equation  discussed  earlier.  To  accomplish  this  we  have  to 
link  the  key  quantities  of  waves  and  radiative  transfer,  viz.,  coherence  function  and 
specific  intensity.  The  relation  between  them  is  obtained  by  computing  the  energy 
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density  using  the  two  concepts.  Thus  we  have 


ie  {(|B„(r)|2>  +  (|Bk(r)|2)}  =  )  f  dSl,I{r,s). 


The  Wigner  transform  provides  us  following  relation: 

1 


(45) 


Defining  X„u  as 


(E^r)E*(r))  =  — j  dk±£^(z,k±). 


1  k'2 

Ip,(z,s)  =  —7-—^cos9£llu(z:kA_), 


(46) 


(47) 


2 rj  (2tt)2 

where  ij  is  the  intrinsic  impedance  of  the  medium,  we  have  from  (45)  and  (46)  /  = 
Xvv  +  Xhh ■  To  facilitate  comparison  with  the  results  of  Ulaby  et  al.  [31],  and  Lam  and 
Ishimaru  [14]  we  will  use  a  modified  version  of  the  Stokes  vector  [12].  Instead  of  the 
standard  form  {/,  Q,  U,  V}  we  use  {(/  +  Q)/2,  (I  —  Q)/ 2,  U,  V}.  Thus,  in  terms  of  X/liy 
defined  in  (47)  our  modified  Stokes  vector  is  {lvv.X)th,  |  (Xvh  +  I hv) ,  —  |  (Ivh  —  Xhv)  }. 

There  is  still  one  difference  that  needs  to  be  ironed  out  before  we  transition  to  the 
RT  equations.  Notice  that  in  our  wave  approach  we  obtained  transport  equations 
for  £,  which  is  the  fluctuating  part  of  the  coherence  function.  On  the  other  hand 
the  phenomenological  RT  equations  are  traditionally  written  for  total  intensities. 
Therefore  we  have  to  express  our  transport  equations  in  terms  of  £.  Notice  that 
£  =  £  +  £,  where  £,  the  average  part  of  £  satisfies: 


([  -  -  III)  }  £“  (z,  kjJ  =  0, 


(48) 


Using  (48)  in  (39)  we  obtain 
-jZ  ~  1  Mk±)  -  ^(kJ-)] }  k^)  = 


\kv 


(2tt) 


da±S>S>%v  ■  ■  v')x 


x$m  <!  kj_  -a±;  i  [gM(kj_)  -  q*( k±)]  -  ^a[qft>{a±)  -  q*u,{ct ±)]  \  £“,u,(z,  a±), 


(49a) 


~JZ  -*Mk±)  -^(k±)]|^(^k±)  =  |^2 


da±S<S<*(v  ■  •  is')  x 
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x$m  jk±  [gM(k±)  -g*(k±)]  -  ^a[g^(a±)  -  g^(a±)]|  £%’(z,a±), 

(49  b) 

Notice  that  this  equation  is  expressed  entirely  in  total  intensity.  Now  we  can  transition 
to  the  phenomenological  RT  equations.  Using  the  relation  between  £  and  /,  we  change 
the  integration  variable  to  solid  angle  and  arrive  at  the  following  equation, 

jcos^+7^1 1?{z,s)=  f  dn'P%a(n}n,)I“(z,s'),  (50  a) 

cosd-^j-  +  j>  if(z,  s)  =  I  dn'pg(n,n’)i?{z,s'),  (so  6) 

where  7  is  the  extinction  matrix  and  P  is  the  phase  matrix.  Implicit  summation 
over  superscripts  a  and  subscript  v  is  assumed  in  (50).  Although  the  structure  of  this 
equation  is  identical  to  that  of  the  RT  (equation  (2)),  the  elements  of  the  phase  matrix 
and  the  extinction  matrices  are  not  the  same  primarily  because  of  coherence  induced 
by  boundaries.  As  mentioned  earlier,  we  assume  that  the  layer  thickness  is  greater 
than  the  mean  free  path  of  the  corresponding  medium.  If,  in  addition,  we  assume  the 
quasi-homogeneous  field  approximation  we  obtain  the  following  expressions  for  the 
extinction  and  phase  matrices: 

7=  cos  6  diag  {2 q",  2 q£,  q"  +  q%,  q"  +  <^}  (51a) 

P$  =  ~k'±-,k'm[acos9  ^  bcosO’]}P“b.  (515) 

The  double  primes  are  used  in  (51a)  to  denote  imaginary  parts.  For  ji  —  v,h, 

75  =  (f*“  ■  v'1)2  75  =  (7  •  h'»)2 

V%  =  (/,“  .  V*)  (M“  •  h»)  7t  =  0  (52) 

Similarly, 

V\fv  =  2  (va  ■  Vb)  (ha  •  Vb )  =  2  (va  ■  h/fe)  (ha  ■  h,b ) 

Vffu  =  (va  ■  Vb)  (ha  ■  tib)  +  (va  •  h/fe)  (ha  •  Vb)  Vtfy  =  0  (53) 
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(54) 


■Py,  =  Kk  =  Pvu  =  0 

Vyy  =  (v“  ■  Vb)  (h“  •  h'&)  -  (Va  •  h'6)  (h“  •  Vb) 

Noting  the  implied  summation  over  a  in  (50)  we  see  that  they  are  identical  to  the 
RT  equations  given  in  Section  2.  Now  we  have  explicit  expressions  for  the  extinction 
matrix  and  phase  matrix  in  terms  of  the  statistical  parameters  of  the  problem  thanks 
to  our  wave  approach. 

We  next  turn  our  attention  to  the  boundary  conditions  (BC).  In  our  wave  approach 
we  obtained  BCs  in  terms  of  ‘stack’  reflection  matrix  R,  whereas  in  the  RT  approach 
the  BCs  are  given  in  terms  of  the  local  interface  reflection  matrices.  We  can  readily 
reconcile  this  apparent  difference.  Note  that  the  BC  in  the  wave  approach  forms 
a  closed  system  whereas  in  the  RT  approach  it  is  ‘open’  (linked  to  adjacent  layer 
intensities).  Let  us  take  a  look  at  the  BC  at  the  (m  —  l)-th  interface.  Rm_2 ,m  can  be 
expressed  in  terms  of  Rm_2,m-i  as  follows, 


R 


m—  l,m  R-m— l,m 


I-  R 


•m— 2,m—  1 


Fm-lR 


•m,m—  1 


-1  ... 


R 


m—2,m—l  r  m—1  m—\,rri' 


(55) 


This  is  the  relation  between  the  stack  reflection  coefficients  of  adjacent  interfaces. 
The  R  and  T  are  local  (single  interface)  reflection  and  transmission  matrices  at  the 
(m  —  l)-th  interface.  On  operating  with  (55)  we  get 


T?d  _  T>  xpw  I  'X'  T?d 

—  m,m—  l^m—  1  * 


(56) 


Notice  that  this  boundary  condition  now  involves  only  local  interface  Fresnel  coeffi¬ 
cients. 

Similarly  we  write  Rm+i,m  in  terms  of  Rm_|_2,m+i  and  hence  obtain  the  BC  at  the 
m-th  interface 


E 


U 

m 


R 


’772+1, m 


E? 


nr 


(57) 


Next  we  take  the  tensor  product  of  (56)  with  its  complex  conjugate  and  average 

w.r.t.  surface  fluctuations.  Employing  the  Wigner  transform  operator,  we  obtain  the 
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boundary  condition  at  the  (m  —  l)-th  interface,  which  is  similar  to  that  of  the  RT 
system.  However,  the  reflection  and  transmission  matrices  used  in  the  RT  system 
correspond  to  the  unperturbed  medium  as  opposed  to  the  average  medium  as  in  the 
case  of  the  statistical  wave  approach.  Similarly  we  can  obtain  the  transport-theoretic 
boundary  conditions  at  the  m-th  interface  using  (57). 

6.  Discussion 

Now  that  we  have  made  the  transition  from  statistical  wave  theory  to  radiative 
transfer  theory  it  is  instructive  to  itemize  the  assumptions  implicitly  involved  in  the 
RT  approach. 

1.  Quasi-stationary  held  approximation. 

2.  Weak  fluctuations. 

3.  Statistical  homogeneity  of  fluctuations. 

These  are  the  three  well-known  conditions  necessary  for  the  unbounded  random 
medium  problem. 

However,  if  the  medium  is  bounded  we  need  to  impose  additional  conditions.  We 
found  that  the  extinction  coefficients  calculated  in  the  wave  approach  and  the  RT 
approach  are  different  and  only  after  applying  further  approximations  can  they  be 
made  to  agree  with  each  other.  The  following  two  additional  conditions  are  required 
for  our  bounded  random  medium  problem: 

4.  Layer  thickness  must  be  of  the  same  order  or  greater  than  the  corresponding  mean 

free  path. 

5.  Weak  surface  correlation  approximation. 

6.  All  fluctuations  of  the  problem  are  statistically  independent. 

In  this  report  we  have  modelled  the  random  media  as  random  continua.  Another 
approach  to  this  problem  is  the  discrete  random  medium  model  [8,  15,  30,  12,  29,  19]. 
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Recently  Mishchenko  et  al.  [19]  (hereafter  referred  to  as  MTL  for  brevity)  derived  the 
vector  radiative  transfer  equation  (VRTE)  for  a  bounded  discrete  random  medium 
using  a  rigorous  microphysical  approach.  This  enabled  them  to  identify  the  following 
assumptions  embedded  in  the  VRTE. 

1.  Scattering  medium  is  illuminated  by  a  plane  wave. 

2.  Each  particle  is  located  in  the  far  field  zone  of  all  other  particles  and  the 
observation  point  is  also  located  in  the  far  field  zones  of  all  the  particles  forming 
the  scattering  medium. 

3.  Neglect  all  scattering  paths  going  through  a  particle  two  or  more  times 
(Twersky  approximation). 

4.  Assume  that  the  scattering  system  is  ergodic  and  averaging  over  time  can  be 
replaced  by  averaging  over  particle  positions  and  states. 

5.  Assume  that  (i)  the  position  and  state  of  each  particle  are  statistically  independent 

of  each  other  and  those  of  all  other  particles  and  (ii)  spatial  distribution  of  the 
particles  throughout  the  medium  is  random  and  statistically  uniform. 

6.  Assume  that  the  scattering  medium  is  convex. 

7.  Assume  that  the  number  of  particles  N  forming  the  scattering  medium  is 
very  large. 

8.  Ignore  all  the  diagrams  with  crossing  connections  in  the  diagrammatic  expansion 
of  the  coherency  dyadic. 

It  is  apparent  that  there  are  distinct  differences  in  the  analyses  for  scattering  from 
discrete  and  continuous  random  media.  Hence  it  is  not  possible  to  make  a  one-to-one 
correspondence  between  the  conditions  of  MTL  and  those  in  this  report.  Below  is  an 
attempt  to  make  a  connection  between  the  two  by  considering  each  condition  derived 
by  MTL  and  relating  it  to  ours.  We  will  denote  the  condition  numbers  derived  by 
MTL  as  MTL  #  and  those  obtained  in  this  report  as  SM 
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MTL  1:-  We  also  have  a  plane  electromagnetic  wave  illuminating  our  system, 
although  as  pointed  out  by  MTL  it  can  be  a  quasi-plane  wave. 

MTL  2:-  We  also  have  implicitly  employed  the  far  field  approximation.  It  is 
embedded  in  SM  1. 

MTL  3:-  This  is  embedded  in  SM  2.  Although  not  explicitly  stated,  the  scattering 
processes  as  mentioned  in  MTL  3  are  avoided. 

MTL  4:-  In  this  report  we  have  restricted  our  attention  to  the  time-independent 
problem  and  hence  did  not  encounter  the  issue  of  ergodicity. 

MTL  5:-  This  condition  is  embedded  in  SM  3. 

MTL  6:-  In  our  problem  we  have  distinct  scattering  boundaries  and  the  character  of 
the  waves  exiting  or  entering  them  are  explicitly  contained  in  the  boundary 
conditions.  Hence  convexity  of  the  scattering  medium  is  not  a  necessary 
condition  for  us. 

MTL  7:-  This  condition  is  embedded  in  SM  4. 

MTL  8:-  This  condition  is  embedded  in  SM  2.  Under  weak  fluctuation 

approximation  we  only  take  into  consideration  the  leading  term  of  the 
intensity  operator. 

Since  the  problem  that  we  considered  in  this  report  involve  scattering  boundaries  we 
have  some  additional  conditions  beyond  those  of  MTL. 

There  are  a  few  more  remarks  that  we  would  like  to  make  before  closing,  (a)  In 
RT  theory  the  medium  is  assumed  to  be  sparse  and  hence  the  “refraction  effects” 
of  the  fluctuations  are  ignored.  Thus  in  the  boundary  conditions  we  should  use 
the  background  medium  parameters  rather  than  the  effective  medium  parameters 
as  derived  in  our  statistical  wave  theory,  (b)  To  arrive  at  (50)  from  (49)  we  have 
ignored  the  contribution  of  evanescent  modes,  (c)  The  condition  about  statistical 
homogeneity  of  fluctuations  may  be  relaxed  by  assuming  it  to  be  statistically  quasi- 
homogeneous  and  we  still  can  arrive  at  our  results  without  much  difficulty,  (d)  The 
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assumption  regarding  underlying  statistics  to  be  Gaussian  is  not  only  a  convenience 
but  also  a  reasonably  good  approximation  in  many  applications.  However,  there  are 
indeed  certain  situations  where  the  statistics  are  not  Gaussian.  Similar  analysis  for 
such  more  general  statistics  are  more  complex  and  involved. 

To  summarize,  we  have  enquired  into  the  assumptions  involved  in  adopting  the 
radiative  transfer  approach  to  scattering  from  layered  random  media  with  rough  in¬ 
terfaces.  To  facilitate  this  enquiry  we  adopted  a  wave  approach  to  this  problem  and 
derived  the  governing  equations  for  the  first  and  second  moments  of  the  wave  fields. 
We  employed  Wigner  transforms  and  transitioned  to  the  system  corresponding  to 
that  of  radiative  transfer  approach.  In  this  process  we  found  that  there  are  more 
conditions  implicitly  involved  in  the  RT  approach  to  this  problem  than  it  is  widely 
believed  to  be  sufficient.  With  the  recent  development  of  fast  and  efficient  algorithms 
for  scattering  computations  and  the  enormous  increase  in  computer  resources  it  is 
now  feasible  to  take  an  entirely  numerical  approach  to  this  problem  without  impos¬ 
ing  any  approximations.  In  spite  of  such  developments,  to  keep  the  size  of  the  problem 
manageable  only  special  cases  have  been  studied  thus  far  [11,  25,  22,  26].  Hence  it  is 
very  much  of  relevance,  interest  and  convenience  to  apply  the  RT  approach  to  these 
problems.  However,  one  should  keep  in  mind  the  assumptions  involved  in  such  an 
approach.  Otherwise  interpretations  of  results  based  on  RT  theory  can  be  misleading. 
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